Structure and phase equilibria of the Widom-Rowlinson model 
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The Widom-Rowlinson model plays an important role in the statistical mechanics of second 
order phase transitions and yet there currently exists no theoretical approach capable of accurately 
predicting both the microscopic structure and phase equilibria. We address this issue using computer 
simulation, density functional theory and integral equation theory. A detailed study of the pair 
correlation functions obtained from computer simulation motivates a closure of the Ornstein-Zernike 
equations which gives a good description of the pair structure and locates the critical point to an 
accuracy of 2%. 

PACS numbers: 61.20.Gy, 61.20.Ja, 64.60.Fr, 05.20.Jj 



I. INTRODUCTION 

The Widom-Rowlinson (WR) model [J is the sim- 
plest model which correctly captures the phenomenology 
of fluid-fluid demixing for systems interacting via short 
range forces and is therefore of fundamental importance 
in the theory of fluids. In particular, the model under- 
goes phase separation at sufficiently high density with 
a critical point which belongs to the Ising universality 
class [1]. The model can be regarded as a non-additive 
hard-sphere mixture in which like species do not interact 
but unlike species exhibit hard-sphere repulsion with a 
given collision diameter a. Although there exist related 
WR-type models in which only the cross interaction is 
non-zero we will reserve the term WR model for that 
with a hard sphere interaction between unlike species. 

Despite the simplicity of the interactions, an accurate 
theory of the bulk structure and thermodynamics of the 
WR model has proved elusive. The lowest order mean- 
field (MF) theory [l|, Q yields a crude description of the 
pair correlation functions and predicts a phase boundary 
between A and .B-rich phases for which the location of 
the critical point is in considerable error when compared 
to recent simulation estimates [1, [|| 0, H, H|. How- 
ever, to go beyond the lowest order theory appears to be 
a very demanding task. The first attempts at system- 
atic improvement were made by incorporating informa- 
tion from higher order virial coefficients into the theory. 
In L 10j virial, activity and cumulant expansions were con- 
sidered. In the virial series was studied for a WR- 
type model in which the AB interaction is that of ori- 
ented hard cubes (chosen to facilitate the calculation of 
higher order virial coefficients). In neither case was re- 
liable improvement obtained. Including additional virial 
coefficients was found to yield either results worse than 
the MF theory or no critical point, due to the appear- 
ance of multiple van der Waals loops in the free energy. 
The origin of these difficulties is the extremely erratic be- 
haviour of the partial sums of the virial series. This issue 
was addressed in [12j in which the virial series of a WR 
type model with Gaussian AB Mayer function was stud- 



ied. Even though the virial series for this model is much 
better behaved than for the WR case, Pade summation of 
terms up to 11th order in density was required to locate 
the critical point with reasonable accuracy. Convergence 
of the virial series in the neighbourhood of the critical 
point for the WR model appears to be much slower than 
for the Gaussian version. Although caution should be 
exercised when drawing conclusions regarding critical be- 
haviour from low density expansions, it seems likely that 
the critical point of the WR model lies outside the radius 
of convergence of the virial series. 

The correlation functions of the MF approximation 
in both the two-component and effective one-component 
versions of the model were studied in Ref. 13. The corre- 
lation functions were shown to display classical behaviour 
at the MF critical point and were proven to be exact in 
the limit of infinitely high dimensionality. A number of 
established integral equation theories have also been ap- 
plied to the WR model [1] although none were able to 
account satisfactorily for the structure or phase bound- 
ary. Of the theories tested only the Percus-Yevick (PY) 
equation displayed a spinodal line, in all other cases (hy- 
perncttcd chain, Martynov-Sarkisov, Rogers- Young and 
Ballone-Pastore-Galli-Gazzillo equations) the theory no 
longer converges in a region of the phase diagram, prior 
to divergence of the partial structure factors. Although 
PY is the best among the standard theories, it still gives 
a poor description of the structure and the critical points 
obtained from the spinodal and virial free energy are 
highly inconsistent. In order to address these difficul- 
ties Yethiraj and Stell (YS) developed an integral equa- 
tion specifically for the WR model in which analytic ex- 
pressions are derived for the direct correlation functions 
between like species caa(t) and css(r) by resumming 
a class of diagrams which can be evaluated analytically 
[14L Il5l | . Unfortunately the YS equation strongly overes- 
timates the structure and the pair correlation functions 
do not compare favourably with computer simulation. 

In Ref. [l6| a fundamental measures density functional 
was constructed for the WR model suitable for applica- 
tion to inhomogeneous situations. The theory predicts 
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TABLE I: Properties of the statepoints at which pair cor- 
relation functions were obtained using computer simulation. 
Shown for each statepoint is the density pA of A particles, 
the fugacity zb of B particles, as well as the corresponding 
density ps of B particles, the composition x, and the total 
density p. 



statepoint 


PA 


zb 


PB 


X 


P 


1 


0.1 


0.8 


0.5815 


0.1467 


0.6815 


2 


0.3 


0.8 


0.3366 


0.4713 


0.6366 


3 


0.4 


0.8 


0.2481 


0.6172 


0.6481 


4 


0.5 


0.8 


0.1758 


0.7399 


0.6758 


5 


0.6 


0.8 


0.1191 


0.8344 


0.7191 


6 


0.7 


0.8 


0.0774 


0.9004 


0.7774 



reasonable behaviour for the pair correlations and im- 
proves slightly on the predictions of the MF theory for 
the phase behaviour. However, the location of the critical 
point remains in error and thus does not permit inves- 
tigation of the structure at statepoints in the vicinity of 
the simulation critical point. Despite all of these efforts 
it can be concluded that the overall level of accuracy 
of the existing theories remains unsatisfactory, given the 
fundamental nature of the model. 

In this work we seek to develop a theory for the WR 
model which provides accurate predictions for the pair 
structure and which is able to locate the critical point to 
an acceptable level of accuracy. In order to base our ap- 
proximations on firm foundations we will use simulation 
results for the pair structure to guide the construction of 
our theory. We focus in particular on the the total cor- 
relation function hij(r) — gij(r) — 1, where gij(r) are the 
radial distribution functions, the bridge function &ij(r) 
and the direct correlation function c.y(r). The paper is 
structured as follows: We begin by discussing the model 
in Section|TTJ In Section HITl we present computer simula- 
tion results for hij (r) , (r) and 6^ (r) , focusing on high 
density statepoints off the critical line where the model 
is least well understood. These simulation results will 
act as motivation for our theoretical approaches in Sec- 
tion IIVI Finally, we summarize our results and suggest 
possibilities for future work. 

II. THE WR MODEL 

The binary WR model is a symmetric mixture consist- 
ing of species A and B. The interaction between like- 
species is ideal, 4>aa{t) = 4>bb(t) = 0, while unlikc- 
species interact via a hard-core potential of diameter 
f 5 ^abW = oo for r < a, and zero otherwise. We 
henceforth take a as the unit of length. Above a cer- 
tain critical density p = (Na + Nb)/V, the WR model 
phase separates into two phases: one phase containing 
predominantly A particles, and the other phase mostly 
B particles. Here, V is the volume, and Na {Nb) 
denotes the number of A (B) particles in the system. 
Due to the symmetry of the model, the compositions 



0.4, 





/ 

/ ■ 




o 




□ 


1 




I.I.I. 



0.2 0.4 0.6 0.8 



X 

FIG. 1: The phase diagram of the WR model in (x,p) rep- 
resentation. The line is the spinodal from the new integral 
equation closure, see Eq. (|12[) . of this work. The circle marks 
the location of the critical point, as obtained in the simula- 
tions of Ref.0. For comparison we also show the critical point 
predicted by the simple mean field theory of [l|, @] (diamond) 
and the density functional theory of Ref. 16; (square). 



of the phases are given by x and 1 — x, respectively, 
with x = Na/(Na + Nb)- At the critical point one 
has x = 1/2. The phase diagram is thus conveniently 
represented in the (x, p)-plane, see Fig. [TJ For densities 
P > Pcrit 7 coexisting phases of composition x and 1 — x 
can be identified. The binodal, which is symmetric about 
the line x = 1/2 in the (x,p) representation, terminates 
at the critical point. In the mean field approximation the 
binodal exhibits a parabolic curvature around the criti- 
cal point (recall the mean-field critical exponent of the 
order parameter (3 = 1/2). The simplest mean-field es- 
timate of the critical density is p C rit = 3/2-7T = 0.4775 
PI 0] , whereas the critical point of the PY spinodal lies 
at /o cr it = 1.12 [I?], [U|]- We emphasize that both of 
these approximations are mean field in character and ex- 
hibit classical critical exponents. In contrast, the best 
current simulation estimates for the critical density are 
p crit = 0.7470(8) @ and p crit = 0.7486±0.0002 [|, which, 
as was pointed out in the Introduction, is not accounted 
for by any theoretical approach. 

A possible source for the discrepancy between simula- 
tion and theory is the fact that the WR model belongs 
to the universality class of the Ising model. For the Ising 
model, (3 ~ 0.326 [HI, implying a flatter binodal. Com- 
puter simulations of the WR model indeed recover Ising 
critical behavior @, However, in order to observe 

the pure Ising exponent /3, the (x, p) representation of the 
binodal is not the most convenient. For Ising systems, 
there is an additional singularity in the specific heat, 
governed by the critical exponent a « 0.109 [l9(. In the 
(x, p) representation, the curvature of the binodal is then 
described by the renormalized exponent /3* = j3/(X — a) 
[20| (in general, critical exponents become renormalized 
if the critical point is approached by varying a quantity 



3 



1.04 




0.76 



0.00 0.15 0.30 0.45 0.60 0.75 0.90 
Pa 

FIG. 2: Phase diagram of the WR model in (pa,zb) repre- 
sentation. The solid curve shows the binodal; the open square 
marks the location of the critical point. The binodal was con- 
structed using simulation data of Ref. |^ combined with finite 
size scaling, and so, on the scale of the above graph, accu- 
rately reflects the true thermodynamic limit form. The closed 
squares (labeled 1 — 6 on the horizontal line in the one-phase 
region), mark the statepoints at which the simulations of this 
work were performed to obtain the pair correlation functions, 
see also Table U 



which is not a field variable). In order to observe the 
pure Ising exponent, the binodal should be represented 
in analogy to the (density,temperature) phase diagram of 
simple fluids. For the WR model, this would be a grand- 
canonical representation, where the density pa = Na/V 
of A species, and the fugacity of B species, are the 
relevant variables 0, l2l| (the choice for A or B is 
of course arbitrary). Shown in Fig. fj] is the phase dia- 
gram in (pa, zb) representation. The curvature of the 
binodal around the critical point, at pa ~ 0.3743 and 
zb ~ 0.93791, is now described by the pure Ising expo- 
nent /3 Q. In contrast to the (x,p) representation, the 
symmetry of the WR model is not obvious from the bin- 
odal of Fig. [2l The symmetry, of course, still exists. In 
the grand-canonical ensemble, it corresponds to the line 
of equal fugacities za — zb- Note that for mean-field 
systems, the curvature of the binodal is not affected by 
the representation, since here a = 0. 



III. COMPUTER SIMULATIONS 

Fantoni et al. have presented computer simulation 
results for gij (r) and Cij (r) for several different values of 
p along the symmetry line x — 1/2. One of the most 
interesting conclusions arising from this work is that the 
Percus-Yevick condition, cab(t) — for all r > 1, is 
satisfied to very high accuracy, even for statepoints ap- 
proaching the critical point. However, it is not at all clear 
whether this property is also maintained off the symme- 
try line; the non-trivial cancellations which apparently 



occur in the diagrammatic expansion of Cy (r > 1) on 
the symmetry line may no longer hold when pa =/= Pb- 
Indeed, a theoretical approach which aims to describe 
phase separation must be able to accurately describe the 
change in pair structure as a function of x. As far as we 
are aware there exists no detailed study of the behaviour 
of the pair structure for x 1/2. 



A. Simulation details 

Motivated by the above considerations, we have per- 
formed simulations for the off-symmetry statepoints 
given in Table Q] To simulate the off-symmetry state- 
points x ^ 1/2, we use a quasi-grand canonical simula- 
tion ensemble, whereby the system volume V, the density 
of A particles pa, and the fugacity zb of B particles are 
fixed, while the number of B particles fluctuates. The 
simulations are performed in cubic simulation boxes of 
edge L = 30, using periodic boundary conditions in all 
d = 3 directions. For the statepoints considered by us, 
this implies approximately 15,000 particles in each simu- 
lation box. To simulate efficiently, a cluster Monte Carlo 
move is used [H, [23|. We specialize to zb — 0.8, which 
is well below its critical value ZB, cr ~ 0.93791 0, ||, 
and inside the one-phase region of the phase diagram, 
see Fig. [2] The density of the A particles is then var- 
ied over the range 0.1 — 0.7. For each statepoint, the 
average concentration ps of B particles is measured, as 
well as the radial distribution functions gij(r). The ra- 
dial distribution functions are evaluated using a standard 
method [24| , and averaged over approximately 5000 inde- 
pendent configurations. For each statepoint, this requires 
an investment of about 120 CPU hours. A total of six 
distinct statepoints is considered. For each statepoint, 
the average concentrations pa and pb, as well as the to- 
tal concentration p, and the composition x, are listed in 
Table H 



B. Analysis of simulation structure 

In Fig. [3] we show the total correlation functions hij(r) 
for three of the considered statepoints. The correlations 
between like species are monotonic and exhibit no sign 
of any oscillatory packing behaviour. The increase of the 
hair) as r — > reflects the tendency of like species to 
overlap in order to maximize the free-volume and hence 
the entropy of the system. Equivalently, this cluster- 
ing behaviour can be viewed as reflecting the attractive 
(many-body) depletion potential acting between spheres 
of species A (£?), induced by the sea of non- interacting 
spheres of species B (A). For example, for small val- 
ues of x we expect hAA(r) — exp[— /3(f>dep(r)} — 1 and 
caa{t) = exp[-/30 dep (r)] +/30 dep (r) - 1, with depletion 
potential 



I3(t>dep{r) 



1 3 1 

1 r H f 

4 16 



(1) 



4 



for r < 2, and zero otherwise. The cross correlation func- 
tion Hab ( r ) is negative for all values of r and indicates 
that in addition to the trivial hard core exclusion there 
is also an effective repulsion between species of opposite 
type. This is a consequence of the clustering of like par- 
ticles and leads to the appealing picture of the bare hard- 
sphere repulsion between species A and B being supple- 
mented by a softly repulsive 'dressed' interaction describ- 
ing particles shrouded by a cluster of like particles. Nat- 
urally, this effective interaction is of statistical origin and 
is therefore not to be taken too literally. The direct corre- 
lation functions are shown in Figs|l]and[5] To obtain the 
Cjj (r) from our simulated hij (r) we apply the method de- 
scribed in Q. The Fourier transform hij(k) yields c\j{k) 
via the Ornstein-Zernike relation. We then construct the 
difference Jij(k) — hij(k) — Cij{k) (the Fourier transform 
of a continuous function in real space) and transform 
back to get Jij(k). We thus obtain a 3 (r) from the differ- 
ence hij(r) — 7<j(r). Both caa{t) and cbb{t) display the 
same monotonic behaviour observed for the total corre- 
lation functions. The ca(r) are shorter range functions 
than the corresponding hu(r), as expected. In Figure [5] 
we show the cross correlation function cab (r) • The form 
of CAs(r) inside the core (r < 1) is quite different from 
the familiar case of an additive hard sphere system for 
which it is found that c a s b(0) < c^ s s (l~) at all densities. 
Outside the core, r > 1, we find that the value of cab{t) 
does not exceed 10~ 3 for any of the simulated statepoints. 
A quantity which is often of interest in liquid state in- 
tegral equation theories are the bridge functions bij{r). 
The bridge functions for a binary mixture interacting via 
pair potentials are defined by the relation 

where cj>ij (r) is the pair potential acting between species 
i and j. In Fig. [6] we show the bridge functions obtained 
from our simulations. Note that data for bij(r) for r < 1 
are not displayed as these are not required for calcula- 
tions of the pair structure for systems with hard core in- 
teractions, see Eq.([2]). While the function bAB(r) is of a 
rather simple form the functions ba(r) take both positive 
and negative values. 

We now consider the implications of the above for con- 
structing approximate theories. The most important in- 
formation to come from the simulations is that the PY 
approximation cab{t > 1) = is satisfied to high ac- 
curacy at all of the considered statepoints. This fully 
supports and extends the findings of Fantoni et al. [8[ 
and implies that the condition cab{t > 1) = should 
be enforced in any approximate theory for this model. 
As the exact core condition Kab^T < 1) = —1 is also 
satisfied, this effectively reduces the problem to that of 
finding accurate approximations for caa{t) and cbb{t)- 
We note that our verification of the PY condition for 
the cross correlation functions invalidates the specula- 
tion in the discussion of Ref. [H. Here, it was suggested 
that a more accurate cab{t > 1) could be achieved by 
modeling the tail with a Yukawa function and adjust- 
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FIG. 3: The correlation functions hij(r) from simulation for 
statepoints 6 (solid line), 4 (broken line) and 2 (dotted line) 
(see Table Hi. Note that h A B{r <!) = -!. 
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FIG. 4: The correlation functions caa(t) and css(r) from 
simulation for statepoints 6 (solid line), 4 (broken line) and 2 
(dotted line). 

ing the free parameters to achieve thermodynamic self 
consistency. This program, while successful for the case 
of hard spheres, would apparently lead to no significant 
improvement over the simple PY approximation for the 
present model. 

The direct correlation functions cu(r) are of signifi- 
cantly shorter range than the ha{r) and possess a rela- 
tively simple monotonic form. It may therefore be useful 
to model the ca(r) by suitably chosen basis functions of 
finite range. By choosing finite range basis functions we 
naturally suppress the development of realistic critical 
exponents (the direct correlation functions are known to 
become long range at the critical point [25]). However, we 
do not necessarily restrict ourselves to mean field critical- 
ity in making this choice. The bridge functions between 
like species ba (r) are very different from those of the hard 
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FIG. 5: The correlation function cab(j) from simulation for 
statepoints 6 (solid line), 4 (broken line) and 2 (dotted line). 
\cab{t > l)j < 10 -3 at all simulated statepoints. 



FIG. 6: The correlation functions bij(r) from simulation for 
statepoints 6 (solid line), 4 (broken line) and 2 (dotted line). 



sphere system [26[ and exhibit regions of both positive 
and negative sign. The function bAB(r > 1), although 
superficially similar to the hard sphere bridge function, 
is found to display a quite different functional form and 
cannot be reasonably fitted using the PY bridge function 
for hard spheres at any effective density. These find- 
ings suggest that modified-hyper-netted-chain (MHNC) 
type approximations [27| , where universality of the hard- 
sphere family of bridge functions is assumed, will not 
prove useful in this case. Indeed, the complex damped 
oscillatory form of the functions bu(r) suggests that ap- 
proaches aiming to directly approximate the bridge func- 
tions should be avoided. 



IV. THEORETICAL APPROACHES 



The simulation results presented in the previous sec- 
tion indicate that a successful theory for the pair struc- 
ture of the WR model may be constructed using the con- 
ditions hAB(r < 1) = — 1, cab(t > 1) = in combi- 
nation with a suitable "ansatz" for the short range func- 
tions Cjj(r). In this section we investigate this possibility. 
Our desire to identify suitable basis functions to describe 
ca (r) leads us to consider a simple virial expansion based 
density functional approximation. This yields analytic 
results for the pair structure which, although only strictly 
valid at low density, actually give a reasonable account 
of the structure over the entire phase diagram. Modifi- 
cation of these results to incorporate the core condition 
leads to a new integral equation closure. 



A. Density functional approach 

Density functional theory (DFT) is a formalism which 
enables the calculation of thermodynamic and structural 
properties of systems subject to spatial inhomogeneity 
|28j |. A key result is the stationarity of the grand po- 
tential with respect to variations in the inhomogeneous 
density fields, 5Q/5pi(r) = 0, where i labels the species. 
Given an explicit functional, this condition yields a set 
of coupled equations for the Pi(v). Ideally, approxima- 
tions within DFT are made directly at the level of the 
free energy, which is a physically intuitive quantity. Cor- 
relation functions of all order can then be generated by 
successive functional differentiation. This is to be con- 
trasted with standard integral equation approaches for 
which the closure is usually introduced at the level of 
the pair correlation functions, and which often make no 
direct reference to, or guarantee the existence of, an ex- 
plicit generating functional. In practice, the distinction 
between integral equation and DFT approaches is fre- 
quently less clear-cut and many approximate DFTs rely 
on correlation functions obtained from integral equation 
theories as input. 

The majority of modern density functional approaches 
are weighted density approximations in which the inho- 
mogeneous density distributions are smoothed by some 
physically motivated set of weight functions [28j ]. Us- 
ing cluster expansion methods [2!| the exact excess 
Helmholtz free energy density can be expressed as a 
power series in the inhomogeneous density fields, pi{r). 
For the WR model, truncation of this series at 0(p 2 ) re- 
covers the ori ginal MF theory [HQ. The only diagram 
contributing to 0(p 3 ) is the triangle diagram consisting 
of two root points and a single field point, all connected 
by Mayer bonds. As the root and field points cannot be 
labelled according to species without either an A-A or 
B-B Mayer bond occuring, the diagram is equal to zero. 



The first correction to the MF theory comes from the 
term 0(p 4 ), which contains only one diagram. By ne- 
glecting terms 0(p 5 ) and higher we obtain the following 
simple approximation 



0T ex [p A ,p B ] 







(3) 



where black and gray field particles are associated with 
the density fields pa(t) and ps(r), respectively, and are 
connected by Mayer bonds. Note that the above dia- 
grams are unlabeled. To convert to labeled diagrams re- 
quires multiplication by the appropriate prefactor (1 and 
1/4). 

At this point we draw on the experience of p revious 
virial expansion studies of the WR model [lj| El, O 03 
and take the truncated expansion ([3]) as the generating 
functional for our correlation functions, at least to a first 
level of approximation. We argue that the above two di- 
agrams contain the dominant structural elements (basis 
functions) for an accurate description of the WR model 
at all densities. Our reasons for this assertion are the fol- 
lowing: (i) investigations of the WR virial series suggest 
that inclusion of higher order diagrams worsens the de- 
scription of the thermodynamic properties, (ii) the pair 
direct correlations and radial distribution functions gen- 
erated from the OZ route give a reasonable account of 
existing simulation results at the simulated statepoints 
(see below), (iii) it can be proven that the MF theory be- 
comes exact in the limit of infinite dimension [§0|. The 
key part of the proof rests on identification of the four 
field particle diagram as the numerically dominant cor- 
rection term to the MF theory. 

It should be emphasized that although Eq. ^ provides 
a reasonable approximation for the thermodynamic func- 
tions over a portion of the phase diagram, the functional 
(Eq.Q) is not a good theory for the phase boundary 
and is not intended as such. Construction of a functional 
which accurately predicts both the bulk binodal and in- 
homogeneous structure is a lofty goal which we do not 
pursue in the present work (see [16| for work in this di- 
rection). Here we present Eq.Q as a means to obtain 
closures at the pair correlation level which may be sub- 
sequently modified and improved. The bulk free energy 
obtained from the uniform density limit of Eq. §3§ is given 

by 



f3F e: 



T^PAPB 



34816 
181440 



^PaPb- 



(4) 



The total Hclmholtz free energy, (3F/V = p\og(p) — p + 
pxlog(x) + p(l — x) log(l — x)+/3F ex /V, displays two van 
der Waals loops as a function of x for sufficiently large 
values of p and is thus unable to account for the demixing 
transition. 

Within DFT the bulk pair correlation functions may 
be obtained using either the test-particle route (minimiz- 
ing the functional in the external field due to a particle 
fixed at the origin) or the Ornstcin-Zcrnike (OZ) route. 
Following the OZ route the inhomogencous pair direct 



Sl.5 

1 

0.5 


-0.5 



v\ AA 


















AB 

iii 


) 


L 2 

r/o 


3 


4 




FIG. 7: The total and direct correlation functions for state- 
point 2. Lines are the results of Eqs.(TBJ)-(r5]). Circles are the 
simulation results. 



correlation functions are obtained by taking two func- 
tional derivatives of the excess free energy functional 



,-(ri,r 2 ) = -/3 



S 2 F ex [{Pi}} 
5pi(ri)6pj(r 2 ) ' 



(5) 



The homogeneous limit is then taken, Cy(ri,r2) — > 
Cij (^12), and the bulk direct correlation functions are sub- 
stituted into the OZ relations to yield the radial distri- 
bution functions. For a binary fluid the OZ relations for 
the homogeneous fluid are given by 



hij(k) 



Cij (k) 



^2pica(k)hij(k), 



(6) 



where the tilde denotes a Fourier transform. Application 
of this prescription to the functional ([3]) generates the fol- 
lowing simple expressions for the bulk direct correlation 
functions 

CA B (r) = f(r)+PAPBf(r)t 2 (r) (7) 
c AA (r) = p 2 B tl(r)/2 (8) 
c BB {r) = p\4(r)/2, (9) 
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FIG. 8: The total and direct correlation functions for state- 
point 2, close to the symmetry line. Lines are the results of 
closure Eq. (fT2)l . Circles are the simulation results. 



where f(r) is the Mayer function, /(r) = —1 for r < 1, 
f(r) = for r > 1, and t\{r) and t^r) are the two lowest 
order chain diagrams given by 



ia(r) 



4/3 1 . . 

r ( 1 -4 r+ i6 r 1 r<2 



(10) 



5040 



(r 3 + 12r 2 + 27r-6)(r-3) 4 



r < 3, 



(11) 



and zero otherwise. Fig. [7] compares the pair correlation 
functions obtained from Eqs.©-© with the simulation 
results for statepoint 2. The level of agreement with the 
simulation gij(r) at such a high density (p w 0.85p C rit) is 
surprising, given the fact that the correlation functions 
are generated from a truncated density expansion. In 
particular, the calculated gAB{r > 1) lies very close to 
the simulation results. Although the general features of 
the Cjj(r) are captured, the overall level of agreement is 
less satisfactory than for the gij (r) . 



B. Imposing the core condition 

The most obvious deficiency of the present approach 
is the violation of the exact core condition, /iab(^) = — 1 
for r < 1. Violation of this condition is a general draw- 
back of the OZ route in DFT studies and only in very 
special cases, e.g. the Rosenfeld functional for additive 
hard sphere mixtures [3lj ] . is the core condition exactly 
satisfied. However, as we are primarily interested in 
the pair correlations this difficulty is easily resolved by 
replacing the closed form expression Eq.(j7j for cab(t) 
with a relation which enforces the core condition. Since 
our simulation results strongly suggest the approxima- 
tion cab { r > 1) = 0, and given that this condition is 
already satisfied by Eq. (J7J , we are led to suggest the fol- 
lowing relations 



h A B{r) 
cab {f) 
caa{t) 
cbb{t) 



= -1 r < 1 
= r > 1 

= PA*i(r)/2. 



(12) 



Combined with the OZ relation, Eq.©, this leads to a 
closed theory for the pair correlation functions. These 
relations correspond to a linearization of the expressions 
for caa{t) and css(r) in the YS integral equation [fH ]. 

In Fig. [5] we show some results obtained using closure 
(fl"2)) at the same statepoint shown in Fig. [71 Eq. (fT2j) was 
solved using standard iterative numerical methods. Im- 
posing the core condition leads to a distinct improvement 
upon the closed form expression Eq.(UJ. The functions 
caa{t) and cbb{t) are identical to those shown in Fig. [Jj 
but cab(t) now li es considerably closer to the simula- 
tion result. The functions qaa{t) and g_Bs(r) remain in 
error for small separations, but the level of agreement 
for r > 1 is improved. There are also small corrections 
to the function Qab{t) over the entire range. We note 
that by imposing the core condition we are effectively 
incorporating many more diagrams (in principle, an infi- 
nite number) into our description of the pair correlations. 
The price we pay for going beyond the simple virial ap- 
proach of §S§ is that we must resort to fully numerical 
solution. 



C. Spinodal line and critical point 

As p is increased for fixed x, the partial structure fac- 
tors Sij(k = 0) diverge at a well defined point. The locus 
of these points defines the spinodal line which divides the 
phase diagram into regions of mechanical stability and in- 
stability. The minimum of this curve, located at x = 1/2 
for the present model, identifies the critical point. It is 
well known that approximate integral equations often fail 
to exhibit a true spinodal but yield instead a no-solutions 
region in the phase diagram within which the theory sim- 
ply fails to converge (see Ref . [32] and references therein) . 
Indeed, the study of Shrew and Yethiraj [J] performed 
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on the symmetry line, and our own investigations for off- 
symmetry compositions, strongly suggest that all stan- 
dard closures, with the exception of the mean-field and 
Percus-Yevick theories, fail to exhibit diverging structure 
factors prior to breakdown of the theory, and are there- 
fore incapable of making any comment regarding the re- 
gion in the vicinity of the critical point. 

In Fig. [TJ we show the spinodal resulting from the clo- 
sure of Eq. (fl2|) . For comparison, we also show the criti- 
cal point from the mean-field theory [J, Q and that from 
the density functional theory of [16j . The critical point 
predicted by Eq. fT!?]) lies remarkably close to the sim- 
ulation result. We find p cr it = 0.762 which compares 
very favourably with the best current simulation esti- 
mates p crit = 0.7470(8) 0] and p cri t = 0.7486 ± 0.0002 
. This represents a substantial improvement upon pre- 
vious theoretical treatments. The numerical solution of 
Eq. (fT2"f for points of high compressibility (i.e. close to 
the spinodal) deserves some additional comment. It is a 
general difficulty of standard numerical methods based 
on Eq.j6|) that, upon approaching the critical point, the 
diverging correlation length renders inadequate methods 
requiring truncation of hij(r) at some finite range R. The 
finite size effects which result from such truncation give 
rise to considerable difficulties when attempting to nu- 
merically assess the critical behaviour of a given integral 
equation [32|. These difficulties have been overcome for 
one-component fluids and mixtures with additive inter- 
actions using specialized algorithms [32|]. However, these 
methods do not generalize easily to non-additive mix- 
tures, such as the WR model, and we have thus resorted 
to more traditional methods of iterative solution [33[ . For 
this reason we make no definite claims regarding the crit- 
ical exponents of the present theory; this would require a 
detailed study using specially tailored algorithms which 
goes beyond the scope of the present work. However, 
the numerical methods we have employed are certainly 
capable of unambiguous determination of the spinodal 
line. This enables us to confirm that the locus of points 
which we have identified is indeed a true spinodal and 
not simply a region of non-convergence. Although we re- 
frain from making final claims regarding the nature of the 
criticality in our equations, we do make the observation 
that the spinodal is distinctly flatter in the vicinity of the 
critic al p oint, compared to the mean- field approaches of 
[E&G3 or the PY theory. This may indicate interesting 
non-classical behaviour, and certainly warrants further 
investigation. 



V. CONCLUSIONS 

Using a combination of computer simulation and the- 
oretical methods we have developed an integral equa- 
tion for the WR model which yields good results for the 
pair structure and predicts the location of the critical 
point to an accuracy of approximately 2%. This rep- 
resents a considerable improvement upon previous the- 
ories which exhibit errors in the range 30 — 50%. Our 
quasi-grand canonical computer simulations provide the 
first detailed information regarding the pair structure 
of the WR model for statepoints off the symmetry line 
(x ^ 1/2) and provide confirmation that the condition 
Cij(r > 1) = is satisfied to a good level of approxima- 
tion over the entire one-phase region. The integral equa- 
tion here developed is very simple to use and requires 
no more numerical effort than solving standard integral 
equations such as PY or HNC. Our choice of basis func- 
tions for Cjj(r) do leave some room for improvement, al- 
beit at the cost of increased numerical effort. A more 
sophisticated scheme could involve basis functions with 
a free parameter, to be determined by enforcing thermo- 
dynamic consistency between virial and fluctuation equa- 
tions of state. Considerable success was achieved in the 
case of additive hard sphere mixtures by constructing an 
approximation for (r) using basis functions taken from 
the low order diagrams in the virial expansion [34] . The 
size of the field particle was treated as a parameter and 
scaled to interpolate between known low and high den- 
sity limits. Whether a similar procedure is also feasible 
for the WR model remains an open question. Following 
completion of this work we were made aware of a very 
recent study in which a triplet level integral equation clo- 
sure was applied to the WR model [35[ . This approach is 
significantly more complicated than that followed in the 
present work but seems to yield very promising results 
worthy of further investigation. 
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